function hydro_init(params)

  % ****************************************************************
  %  output info
  % ****************************************************************

  % summary
  disp(sprintf('Grid Size   : (%d, %d)', params.Nx, params.Ny));
  disp(sprintf('Grid Length : %f', params.Lx));
  disp(sprintf('Cell Length : %f', params.Dx));
  %disp(sprintf('Time steps  : %d',  params.t));
  %disp(sprintf('Hydro Size  : %6.2fMB', 3*params.Nx.*params.Ny.*5.*params.Tsteps.*64/8/1024/1024));



  % ****************************************************************
  %  prepare plotting
  % ****************************************************************

  if (params.plotHydro)
    % prepare figures for plotting
    close all;

    % load color maps
    load('cm_turb', 'cm_turb');
    cm = cm_turb;

    global cmSize;
    cmSize = 512;

    % figures
    fig  = figure; 
    set(fig, 'WindowStyle', 'docked');
    set(fig, 'ColorMap', cm);
    set(gcf, 'renderer', 'zbuffer');

    refresh();
    pause(0.001);
  end
  

  % ****************************************************************
  %  compute spectral derivative kernels
  % ****************************************************************

  tic;

  fprintf('computing spectral derivative matrices...  ');

  Nx = params.Nx;
  Ny = params.Ny;
  Dx = params.Dx;
  Dy = params.Dy;

  % can upgrade on this -- trefethen pg 6

  % the derivative code assumes Nx = Ny for now. 
  % Looks like that might in fact stay this way forever.

  global d1x d1y d2x d2y;
  global gd1x gd1y gd2x gd2y;

  if size(d1x) ~= [Nx, Nx] | size(d1y) ~= [Ny, Ny]
    d1x = zeros(Nx, Nx);
    d1y = zeros(Ny, Ny);
    d2x = zeros(Nx, Nx);
    d2y = zeros(Ny, Ny);

    % x derivatives
    for x = 1:Nx
      for y=1:Nx
        s1 = 0;
        s2 = 0;
        for k = -Nx/2 : Nx/2
          s1 = s1 + 2. * pi * i / Nx / Nx * k * exp(2 * pi * i * k * (x - y) / Nx);
          %s2 = s2 + (2. * pi * i / Nx * k)^2 / Nx * exp(2 * pi * i * k * (x - y) / Nx);
        end
        d1x(x,y) = real(s1)/Dx;
        %d2x(x,y) = real(s2)/Dx/Dx;
      end
    end

    % y derivatives
    for x = 1:Ny
      for y=1:Ny
        s1 = 0;
        s2 = 0;
        for k = -Ny/2 : Ny/2
          s1 = s1 + 2. * pi * i / Ny / Ny * k * exp(2 * pi * i * k * (x - y) / Ny);
          %s2 = s2 + (2. * pi * i / Ny * k)^2 / Ny * exp(2 * pi * i * k * (x - y) / Ny);
        end
        d1y(x,y) = real(s1)/Dy;
        %d2y(x,y) = real(s2)/Dy/Dy;
      end
    end

  end

  tder = toc;
  fprintf('done in %.2f sec\n', toc);
end
